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Abstract. Using molecular dynamics simulations, we study the slow dynamics of 
supercooled liquids confined in a random matrix of immobile obstacles. We study 
the dynamical crossover from glass-like to Lorentz-gas-like behavior in terms of the 
density correlation function, the mean square displacement, the nonlinear dynamic 
susceptibility, the non-Gaussian parameter, and the fragility. Cooperative and spatially 
heterogeneous dynamics are suppressed as the obstacle density increases, which lead 
to the more Arrhcnius-like behavior in the temperature dependence of the relaxation 
time. Our findings are qualitatively consistent with the results of recent experimental 
and numerical studies for various classes of spatially heterogeneous systems. We also 
investigate the dependence of the dynamics of mobile particles on the protocol to 
generate the random matrix. A reentrant transition from the arrested phase to the 
liquid phase as the mobile particle density increases is observed for a class of protocols. 
This reentrance is explained in terms of the distribution of the volume of the voids 
that are available to the mobile particles. 
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1. Introduction 

The transport properties of fluids in a heterogeneous environment are of great 
importance in physics, chemistry, and engineering [H |2]. These systems include fluids 
confined in walls, standing thin-film liquids, and fluids adsorbed in random media. 
The effect of spatial confinement is especially important to the study of the glass 
transition of supercooled liquids. Recent experiments and computer simulations have 
revealed that the glass transition temperature, transport properties, and microscopic 
dynamics sensitively change in the presence of spatial confinement [3]. Moreover, it 
is expected that an understanding of these phenomena may lead to a deeper insight 
into the growing length scales of the cooperative motion of atoms, which escorts the 
drastic slowing down of dynamics near the glass transition point [HE]. A fluid in 
randomly distributed immobile obstacles is an ideal model to study the effects of 
the confinement on the glassy slow dynamics, in fact, a fluid in random media is 
interesting in its own right. This system is introduced as a model system of a crowded 
environment, and the transport phenomena contained within has attracted a great deal 
of attention in the biophysics community [HI El [HI [9] . This system can also be seen as 
a generalization of the Lorentz gas problem to the multi-particle case pUl EED, H21 EE]. 
Furthermore, this system can be regarded as a model of binary mixtures with a disparate 
size ratio [HI HHJ [161 HZl HEl UHl EQl EH [22l [23l El]- The immobile obstacles in the 
random media are interpreted as large particles in the binary liquids because of the 
huge asymmetry in time scales between the small and large particle components. 

In recent years, various molecular dynamics (MD) simulations have been performed 
to study the slow dynamics of the fluids in random media (251 12H1 [271 123 [291 Ell E] . 
Theoretically, the slow dynamics of mobile hard spheres in the presence of the immobile 
hard spheres of the same size has been intensively investigated using the replica method 
combined with the mode-coupling theory (RMCT) [321 [331 EH ESj- These studies have 
examined the dynamic phase transition from liquid to non-ergodic arrested states. Two 
notable results were predicted. The first is that the slow dynamics can be characterized 
by two types of dynamics: Type A and Type B dynamics [36J . When the mobile particle 
density, p m , is large and the immobile particle density, p i; is small, the system undergoes 
a conventional glass transition, in which the onset of slow dynamics is signaled by the 
discontinuous emergence of two-step relaxation in the density correlation function. This 
is referred to as Type B transition. As the immobile density pi increases, the glass 
transition point of the mobile particles decreases drastically. At an even larger pi, the 
dynamics qualitatively changes; one-step slow relaxation sets in at large wavelengths, in 
which the tail of the relaxation curve continuously grows and progressively propagates 
toward the smaller wavelengths as pi increases, until the dynamics of the mobile particle 
freezes. This freezing is a localization transition due to blocking by the percolating 
network of the immobile particles. This is called Type A transition. The crossover 
from Type B to Type A transition is reminiscent of the behaviors observed in many 
heterogeneous systems such as the glass-to-gel crossover of attractive colloids [371 [38] an d 
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Figure 1. Schematic illustrations of two protocols to generate the random matrix: 
(a) Quenched-Annealed (QA) and (b) Equilibrated Mixture (EM) protocols. 



the glassy slow dynamics of binary mixture systems with a disparate size ratio [TBI EE] • 
The second prediction is the existence of a reentrant transition in the small p m regime, 
in which the arrested mobile particles melt as the mobile particle density p m increases 
at a fixed pi. 

To verify these theoretical predictions, we [391 SO] and Kurzidim et al. [TT1 IT2] have 
independently and numerically investigated fluids in immobile obstacles. All of these 
studies have confirmed that there is a crossover from Type B to Type A transition as 
the immobile particle density increases. These studies also found no reentrant transition 
for the system that was studied by RMCT. 

We have also found that the dynamic arrest line (or more precisely iso-relaxat ion- 
time line) sensitively depends on the protocol that is used to generate the configuration 
of the randomly distributed immobile particles [391 SO] • Two types of these protocols 
have been studied; the first is the Quenched-Annealed (QA) system, in which the 
immobile particles are initially equilibrated before their positions are quenched (see 
Fig. DJa)). Mobile particles are then inserted into the system and their dynamics are 
monitored after equilibration. This QA protocol offers a natural choice to model the 
experimental setups of random media (such as porous materials) and has been adopted 
in the theoretical analysis of RMCT [32} |33l EH [35] and in numerical studies by Kurzidim 
et al. [UJ H2] • The second protocol studied is the Equilibrated Mixture (EM) protocol, 
in which all of the particles are run in the simulation box and, after their equilibration, 
the motions of a fraction of the particles are quenched (see Fig. H^b)). The dynamics of 
the mobile component is monitored after waiting long enough for the mobile particles 
to equilibrate in the presence of the immobile particles. This protocol is appropriate 
as a model of the dynamics of the fast (small) particle component in binary mixtures 
with a disparate size ratio, in which the two time scales of each component are well 
separated. We found no reentrance for the system prepared with the QA protocol, but 
surprisingly, we did observe reentrance using the EM protocol [391 HQ], ft should be noted 
that this reentrance is distinct from that predicted by RMCT. It was speculated that 
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this reentrance can be attributed solely to the configurations of the immobile particles 
prepared by the EM protocol; the configuration of immobile particles is automatically 
"optimized" to provide more pathways for the mobile particles than those prepared in 
the absence of mobile particles (using the QA protocol). 

In the present paper, we investigate the dynamical properties of the fluids in random 
obstacles (which were also examined in our previous papers [39], HQ]) in more detail, 
focusing on various quantities that characterize the slow dynamics near the glass and 
the localization transition. We evaluate the nonlinear dynamic susceptibility, the non- 
Gaussian parameter, and the fragility across the entire parameter space of (pi,p m ). We 
also quantify how the dynamics of mobile particles sensitively depends on the protocol 
used to generate the random matrices. 

This paper is organized as follows. In Sec. [2j we briefly review our model and 
simulation method. In Sec. EJ the numerical results are given; we first describe 
how the dynamics changes from Type B to Type A by calculating various dynamic 
quantities. In the latter subsection of Sec. [3] we discuss the sensitivity of the dynamics 
of mobile particles to a geometry of random configurations of the immobile particles by 
calculating the pore-size distribution. In Sec. EJ we summarize our results and provide 
our concluding remarks. 

2. Model and methods 

We perform MD simulations for two types of systems: a binary mixture interacting 
with the soft-core potential and a monodisperse hard sphere. The binary mixture is 
used to explore the entire parameter space of (pi, p m )] bidispersity is required to avoid 
crystallization at the small obstacle density limit — > 0. The monodisperse hard sphere 
is used to investigate the region at which the crossover from Type B to Type A transition 
takes place. This is also the region where the reentrant transition was predicted by 
RMCT. Monodispersity is of importance to explore the physical mechanism near the 
crossover without the risk of it being obscured by the softness of the potential or by the 
bidispersity of the system. 

The binary soft-core mixture consists of an equal number of two types of particles 
with the total number N = N\ + N 2 = 500 + 500. They interact via the soft-core 
potential 



where a a /3 = (o~ a + up)/2 and a, (3 6 {1, 2}. The size and mass ratio were 02/01 = 1.2 
and m 2 /mi = 2, respectively. The total number density was fixed at p = (Ni+N 2 )/ L 3 = 
0.8crf 3 , in which the system length was L = 10.7701 under periodic boundary conditions 
(PBC). The units of length, time, and temperature were considered to be 01, \jm\o\je, 
and e/ks, respectively. For each simulation run, Ni particles were picked up from 
the N particles randomly and fixed their positions. N m = N — Ni particles were left 
mobile. In this model, we used the number densities defined by pi = Nip cS /N and 
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Figure 2. Dynamic phase diagram of the binary soft-core mixture generated by the 
EM protocol, pi is the immobile (obstacle) particle density, and p m is the mobile 
(fluid) particle density. The arrested states are defined as the points beyond which the 
a-rclaxation time r a exceeds 10 3 . 



p m = N m p e g/N as the system parameters, in which p cff = N{e/kBT) l / A a c ^ /V is the 
effective density of this model [13]. Here, a e fj is the effective particle diameter defined 
by 0cff 3 = J2a,p=i,2 x a%i30~ai3 3 , where X\ = Ni/N = 1/2 and x 2 = N 2 /N = 1/2 are the 
mixture compositions. The states that were investigated here were as follows: p c s = 0.5, 
0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.15, 1.2, 1.3, 1.4, and 1.45. The corresponding temperatures 
were as follows: T = 21.61, 10.42, 5.624, 3.297, 2.058, 1.350, 0.992, 0.772, 0.651, 0.473, 
0.352, and 0.306, respectively. We controlled pi and p m by changing iVj (or N m ) and 
p c fj (or T). The number of mobile particles was chosen to be N m = 10 ~ 900. The 
velocity Verlet algorithm was used to integrate Newton equations with time steps of 
0.001 ~ 0.005. 

The monodisperse hard sphere system includes N identical hard spheres with mass 
m and diameter a in a cubic box of volume V under PBC a and ^ma 2 /ksT were used 
as the units of length and time, respectively. The temperature was fixed as ksT = 1. 
The standard event-driven algorithm was used for particle collisions [S]. The number 
densities pi = Ni/V and p m = N m /V were controlled by changing N iy N m , and V. 
For both systems, we carefully checked the system size dependence and the sample 
dependence of the observables throughout the study. Two types of protocols, the QA 
and EM protocols, were employed to generate the random matrices. 



3. Numerical results 

3.1. Dynamic phase diagram 

We first determined the dynamic phase diagram for the whole (p,, p m )-space by 
performing MD simulations of the EM systems of the binary soft-core mixture. In 
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Figure 3. Various timc-depcndcnt quantities for the binary soft-core mixture 
generated by the EM protocol at the four state points (A) (Ni,p e g) = (100,1.45), 
(B) (Ni,p eS ) = (300,1.3), (C) (N hPeS ) = (500,1.15), and (D) (N hPeS ) = (900,0.7), 
as denoted in Fig. [2] (a) the self-part of the intermediate scattering function F s (k,t) 
with k = 2ir, (b) the mean square displacement (Ar 2 (i)), (c) the nonlinear dynamic 
susceptibility X4(/c,i) with k = 2tt, and (d) the non-Gaussian parameter ct2(t). 



Fig. |2l the dynamic phase diagram is plotted as a function of pi and p m . The dynamic 
arrest line is defined as the points at which the a- relaxation time r a reaches 10 3 . We 
confirmed that varying the criteria for r a simply shifts the dynamic arrest line back and 
forth, but that its qualitative behavior remains intact. r a is determined by calculating 
the self-part of the intermediate scattering function for mobile particles, 

1 l Nm - \ 
F.(M) = jj~ (Y,eMik ■ {r 3 {t) - r y (0))]^ , (2) 

where k is the wave vector, k = \k\, and rj(t) is the position of the j-th particle. We 
defined r a by F s {k = 2ix, r Q ) = 0.1. 

As pi increases, the dynamic arrest line (or the glass transition points) drastically 
decreases. These features were well documented in the previous simulations studies [281 
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M> EDI EI] and in RMCT [321 E31 EH This tendency is sustained up to the region, 
beyond which a small reentrant pocket is observed. This reentrance will be discussed in 
a later section. 

3.2. Intermediate scattering function 

The numerical results of conventional dynamical quantities displayed in Fig. [3] 
demonstrate how the immobile particle density pi affects the relaxation processes of 
mobile particles. In Fig. El^a), the time evolution of F s (k,t) with k = 2tt is plotted for 
the four state points (A)-(D) that are indicated in Fig. [2j at which the a-relaxation 
times are almost the same. The figure clearly indicates that there are two types of 
distinct dynamics depending on pi. 

In the small pi (immobile particle density) regime pi <C 0.5, F s (k,t) exhibits two- 
step relaxation with a well-developed plateau, which is a hallmark of slow dynamics near 
the glass transition point. We found that the shoulder of the plateau discontinuously 
appears as one approaches from the fluid side to the arrested phase [391 ED]- This 
behavior is typical of slow dynamics near the glass transition point and is referred 
to as Type B transition in the MCT community. However, as the mobile particle 
density p m decreases and pi increases, the relaxation profile of F s (k,t) becomes quite 
different from that of Type B transition, i.e., F s (k,t) shows a single step relaxation 
with a long tail (see F s (k,t) at the state (D)). It is also observed that the amplitude of 
the tail increases continuously as one crosses the arrested phase and that this increase 
incipiently starts from the lowest wavelength and propagates to the shorter scales as 
Pi increases [39], HQ]. This behavior is known as the hallmark of Type A transition (or 
the localization transition) as predicted by RMCT [321 E31 EH ES] and demonstrated 
by simulations for various spatially heterogeneous systems, such as binary mixtures of 
large and small particles [HI [17] and colloidal gels [37J 138] . 

3.3. Mean square displacement 

The qualitative change from Type B to Type A dynamics is also observed in the mean 
square displacement (MSD) for mobile particles, 



The results are plotted in Fig. EJb). It is known that in the Type B dynamics, the MSD 
exhibits a plateau at the /3-relaxation time regime where the plateaus are observed for 
F s (k,t), in which the tagged particle is trapped by surrounding particles. However, at 
the small p m limit (see (D) of Fig. EHb)), anomalous subdiffusive behavior Ar 2 (t) ~ t a 
(a < 1) is observed, in which the system is almost Lorentz-gas-like. On a short 
timeframe, the mobile particles at the state point (D) can explore longer distance than 
those at (A)-(C), because the mobile particle density is small. At a longer timeframe, 
however, the diffusion becomes very slow due to the developing network of immobile 
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particles that hinders the ability of the mobile particles to explore the long distance; 
this system can be explained in terms of the percolation theory [15]. The subdiffusion 
exponent a ~ 0.3 is consistent with 0.32 predicted for the Lorentz gas [TQl HH C21 fl3j . 

3.4- Nonlinear dynamic susceptibility 

We next investigate the nonlinear dynamic susceptibility or four-point correlation 
function for mobile particles, X4{k,t), which is a measure that is used to quantify the 
extent of the dynamic heterogeneities H7J HSj . X4,{k,t) is defined as the variance of 
the fluctuations of the self part of the intermediate scattering function by 

X i(k,t) = N m [{P?(k,t)) - (F s (k,t)) 2 ]. (4) 

Here F„(k,t) = (F s (k,t)) and 

In the literatures (17], iV" 1 X)i=i cos(k ■ Ar^(t)) is conventionally used as the definition 
of F s (k,t). Under this definition, x±{k,t) decays to a constant 1/2 at t — > oo. However, 
as we demonstrate here, the peak of XA{k,t) for the Type A regime grows more mildly 
than it does for the Type B regime. To demonstrate this suppression of the peak 
and thus the dynamic heterogeneities in the Type A regime without the results being 
obscured by a constant plateau at a large t, we have adopted an alternative definition of 
F s (k, t) by taking the average over the angular components of the wave vector k, which 
leads to eq.(|H]). Note that both definitions of F s (k,t) lead to an identical averaged 
value F s (k,t) = (F s (k,t)) due the isotropic nature of the system, but that the new 
definition removes the unwanted constant for x±{k,t) at t — > oo. In Fig. [3]^c), the time 
evolutions of x<±{k, t) are plotted for the four state points (A)-(D). In the Type B regime, 
X4,{k,t) exhibits behavior that is typical for bulk glass, i.e., a pronounced peak at the 
a-relaxation time, whose height grows rapidly as the density increases and is preceded 
by algebraic growth in the /3-relaxation regime. In the Type A regime at state (D), 
however, Xi{k,t) neither grows nor shows a strong peak, even after a long period of 
time. This results implies that dynamic heterogeneities play a minor role in the slow 
dynamics of this regime. Similar behavior of X4,{k, t) has been reported for colloidal gels 
whose slow dynamics are caused by geometrical constraints [491 150] . 

3.5. N on- Gaussian parameter 

The non-Gaussian parameter (NGP) is another typical quantity that is suitable to 
monitor the effect of the heterogeneities inherent in the system. We calculated the 
NGP a 2 (t) defined by 

, , 3(Ar 4 (t)) 
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Figure 4. (a) The a-relaxation time, r a , as a function of the inverse temperature 
T g /T, for several A^'s, N, = 0, 100, 200, 500 and 800. T g is defined as the temperature 
at which r a = 10 4 . (b) ^-dependence of the fragility index m. 

where (Ar 4 (£)) = (l/N m )(J2f=i\rj(t) — r*,(0)| 4 ). a>2(t) reveals how the distribution 
of the single-particle displacement, \Afj(t)\, at time t deviates from the Gaussian 
distribution [51]. The profiles of are indicated in Fig. [3]^d). 

It is observed that cti (t) develops and exhibits the pronounced peak before the a- 
relaxation time in the Type B regime. It is known that an increase in the maximum 
a2(t) is synchronized with the growing dynamic heterogeneities near the glass transition 
point [52]. As the immobile density p^ increases, the height of the peak decreases. At 
the largest pi, the point (D), it is hard to see the peak. This trend is qualitatively 
similar to that of the nonlinear dynamic susceptibility, Xi(k>t)- Recently, Flenner and 
Szamel have proposed a new non-Gaussian parameter (NNGP) [53]. They argued that 
the conventional NGP is more susceptible to particles moving faster than those moving 
more slowly, whereas NNGP is more susceptible to slowly moving particles. It would be 
beneficial to compute and compare the NGP, NNGP, and x^ik, t) in order to corroborate 
the role of the dynamic and static heterogeneities in the confined systems. 



3.6. Fragility 

The "fragility" is a concept to quantify the deviation of the temperature dependence 
of the viscosity, diffusion coefficient, and the relaxation time from the Arrhenius 
behavior The fragility index m is commonly used as a measure of the fragility 

and is defined by the steepness of the increase of r a upon decreasing the temperature; 

d{T g /T) 

where T g is the glass transition point, m depends on the material of the glass formers. 
Generally, fragile liquids with large m's tend to exhibit more pronounced and more 
temperature-sensitive dynamic heterogeneities than do more Arrhenius-like fluids with 



m 



(7) 
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smaller m's (or stronger liquids). As was demonstrated in the previous sections, the 
dynamic heterogeneities are suppressed as the immobile particle density pi increases. 
Therefore, it is natural to expect that the system concomitantly becomes stronger (more 
Arrhenius-like). We examined the temperature (or the effective density p c s) dependence 
of Tc's for various N^s (the number of the immobile particles) for a binary soft-core 
mixture. Note that we have used p c s to control T. The change of p c ^ changes pi (and 
also p m ) slightly, and this shift may make it difficult to quantify the effects of the fixed 
number of obstacles on the fragility. We believe, however, that this effect is negligible 
for the range of the temperatures that we have explored. 

The "Angell-plot" , the log 10 r Q -vs-l/T plot, of our system is shown in Fig. 11(a), 
in which the temperature T is scaled by the "glass transition temperature" T g . T g 
has been defined by the point at which r Q reaches 10 4 that is slightly longer than the 
criteria used to draw Fig. [2j The dependence of the fragility index m on the number of 
immobile particles iVj is plotted in Fig. H^b). We observe that m is remarkably sensitive 
to the density of immobile particles. The fragility index is m ~ 14 for bulk glass but it 
decreases to m ~ 4 at the largest density of immobile particles. Our observation differs 
from the one reported for a binary mixture system with large size ratios (< 3) by Kurita 
et al. [23]. They indicated that the fragility changed non-monotonically, though mildly, 
by changing the size ratio of small and large particles and their densities. It would be 
interesting to study how this trend changes as the size ratio of the two components 
increases, resulting in the time scales for each component becoming decoupled. 

We remark that the fragility, which behaves in a manner that is qualitatively 
similarly to ours, has been experimentally obtained in polymeric systems confined in 
porous media recently [55], in which the crossover from a non-Arrhenius to Arrhenius 
temperature-dependence of the relaxation time was observed as the pore size became 
smaller. It was speculated that the decrease of the fragility as the effect of the 
confinement is enhanced is universal and should be observed for other types of confined 
systems such as those with the solid-liquid or air-liquid interfaces [56] . 

Finally, it should be noted that even the largest values of m reported here are 
still very small compared with conventional molecular systems [57J. This discrepancy 
occurs because the glass transition temperature T g defined above is far higher than those 
observed for real glasses, which is due to the limited time windows that the simulations 
can access. It is noteworthy that our simulation results still exhibited qualitatively 
similar behavior for m as the experimental results, despite the large time-scale differences 
between them. 

3. 7. Reentrant transition 

In this subsection, we examine the effect of the configurations of immobile particles on 
dynamics of the mobile particles by comparing results obtained with the EM protocol 
to those obtained with the QA protocol in the monodisperse hard sphere system. We 
used the one-component hard spheres to study the mechanism behind the configuration- 
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Figure 5. Dynamic phase diagram of the monodisperse hard sphere system with a 
random matrix generated by the two different protocols: (a) QA and (b) EM. 




Figure 6. The p m dependence of the mean square displacement, (Ar 2 (t)), at the 
immobile particle density (a) pi = 0.3, (b) 0.43, and (c) 0.5 of the QA systems 




Figure 7. The p m dependence of the mean square displacement, (Ar 2 (t)), at the 
immobile particle density (a) pi = 0.3,(b) 0.43, and (c) 0.5 of the EM systems. 



dependence of dynamics near the Type B- Type A crossover without the results being 
obscured by the softness of the potential or by the bidispersity of the system. We 
calculated the MSD of the mobile particles and determined the dynamic phase diagrams 
for both the QA and EM systems. The results are plotted in Fig. [51 Here, the dynamic 
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Figure 8. The pore-size distribution P(v poIe ) for various mobile particle densities 
p m at the fixed immobile particle density pi = 0.5 of the EM systems. 

arrest line is determined as the points at which the MSD reaches 10 2 in the simulation 
time t = 10 4 . Figs. [6] and [7] show the dependence of the MSD on p m at several p^s for 
both the QA and EM systems. 

As indicated by Fig. [5(a), no reentrant transition is observed for the QA system. 
The dynamic arrest line monotonically decreases as pi increases, which is compatible 
with the recent numerical simulations for the related QA systems jH}H2]. Indeed, Fig. [6] 
indicates that the slope of the MSD monotonically decreases as p m increases at a fixed 
Pi. On the one hand, this result is hardly surprising because the slowing down of the 
mobile particle dynamics is mainly due to the geometrical confinement by the immobile 
particles. On the other hand, the EM system shows the reentrant pocket at a finite 
p m , which is clearly seen in Fig [5(b). The similar reentrance pocket has been observed 
in the binary soft-core mixture (see Fig. [2]). The dynamics of the mobile particles are 
accelerated in spite of the increase in p m . This reentrance can be clearly seen in Fig. [7(b) 
and (c); by increasing p m at the fixed large p i: the slope of the MSD at long times first 
increases and then gradually decreases. 

In our previous study [391 HO], we have speculated that the origin of this reentrance is 
due to the change of the equilibrium structure of the immobile particles in the presence 
of the mobile particles, which are equilibrated together when the random matrix is 
generated. To verify this speculation, we investigated the distribution of the pore size (or 
the free volume available for the mobile particles) generated by the immobile particles. 
The pore-size distribution is determined as follows. Using a three dimensional Delaunay 
triangulation algorithm, the total space of the system is divided into non-overlapping 
tetrahedrons. The vertices consist of the positions of the immobile particles. The volume 
distribution of the tetrahedrons, P(v poTe ), for EM systems is computed. If the volume 
of the tetrahedron is much larger than that of the particle, v porc > 7ia 3 /6 ~ 0.52a 3 , the 
mobile particle can access the pore. The available pore-sizes for the mobile particles are 
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not identical to the available pathways that are available to them, but their distribution 
function still provides reliable information on the dynamics of the mobile particles in 
geometrical confinement. As observed in Fig. EJ the height of the tail of P(v porc ) at 
Vpore > 0.52cr 3 increases as p m increases at a fixed value of pi = 0.5 monotonically. 
This result indicates that the free volumes available for the mobile particles increases, 
which leads to the reentrant behavior of the MSD that is observed in Fig. [7(c) . Note 
that the tails monotonically increase (at least up to p m = 0.4), but the dynamics slow 
down again around p rn ~ 0.2, as indicated by Fig. [5(b). This result occurs because the 
glassy dynamics of the mobile particles sets in while the localization effect (due to the 
geometrical confinements) becomes smaller. 

These results quantify our speculation that the immobile particles adjust themselves 
during an equilibration run to prepare for the presence of the mobile particles so that 
the free volumes for both components are entropically maximized and leave more 
available geometrical spaces (and pathways) for the mobile particles, delaying the 
percolation transition to larger values and thus leading to faster dynamics of the mobile 
particles. The sensitivity of the percolation point to the protocols used to generate the 
matrix configuration has been previously studied in several contexts and our results are 
consistent with these results [291 USB I3"T] . 

We also speculate that counterintuitive effects similar to the reentrance discussed 
above are prevalent in systems in which the disparate time scales are entangled. A 
binary mixture consisting of large and small particles studied by MD simulations may 
be a good example. Voigtmann et al. have numerically studied the dynamics of such a 
binary mixture and have shown that the diffusion of small particles becomes slower when 
the interactions between small particles are turned off j22j [23]. One may speculate that 
this observation is relevant to our finding of the reentrant pocket for the EM system. The 
turn-off of the interactions makes small particles "invisible" to each other and makes 
the large particles behave as if there are fewer small particles around them, which 
makes the configuration of the large particles become more QA-like than those with full 
interactions. This effect might be difficult to observe via standard static quantities (such 
as the static structure factors) but may be detected easily via the pore-size distribution 
function. An accurate theoretical method to evaluate the subtle differences in the static 
structure factors would be desirable to investigate how the protocol dependence changes 
the dynamic phase diagram using RMCT [58] . 

4. Conclusions 

In this paper, MD simulations have been performed to examine the dynamical properties 
near the arrest points of simple fluids confined in random media. We calculated various 
quantities for the whole range of the mobile and immobile densities, including the 
intermediate scattering function, the mean square displacement, the nonlinear dynamic 
susceptibility, the non-Gaussian parameter, and the fragility. We found that all of these 
quantities exhibited qualitative changes as the density of the mobile/immobile particles 
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was varied. At the limit of small immobile particle density, all of the observed quantities 
indicated the strong signs of dynamic heterogeneities near the glass transition point, such 
as the enhanced nonlinear dynamic susceptibility, the increased peak heights of the non- 
Gaussian parameter, and the fragile behavior of the relaxation time in its temperature 
dependence. At the opposite limit, in which a small number of mobile particles diffuse 
amid abundant obstacle particles, the signs of the dynamic heterogeneities were all 
suppressed. The nonlinear dynamic susceptibility and the non-Gaussian parameter 
exhibited no peak, and the temperature dependence of the relaxation time was almost 
Arrhenius. To understand the underling physics behind the reentrant transition near the 
Type B to Type A crossover, we have carefully quantified how the statistical properties 
of configuration of the random matrix can be altered by the different protocols used 
to generate them by calculating the pore-size distribution generated by the immobile 
particles. Throughout this paper, we refer the change from Type B to Type A dynamics 
as the crossover. According to RMCT [321 [331 EU [35] , this change should be associated 
with a higher order MCT transition. However, this transition is too subtle to be clearly 
observed at the resolution of the current simulations. 
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